Robust Recovery of Optimally Smoothed Polymer Relaxation Spectrum from Stress Relaxation Test Measurements

The relaxation spectrum is a fundamental viscoelastic characteristic from which other material functions used to describe the rheological properties of polymers can be determined. The spectrum is recovered from relaxation stress or oscillatory shear data. Since the problem of the relaxation spectrum identification is ill-posed, in the known methods, different mechanisms are built-in to obtain a smooth enough and noise-robust relaxation spectrum model. The regularization of the original problem and/or limit of the set of admissible solutions are the most commonly used remedies. Here, the problem of determining an optimally smoothed continuous relaxation time spectrum is directly stated and solved for the first time, assuming that discrete-time noise-corrupted measurements of a relaxation modulus obtained in the stress relaxation experiment are available for identification. The relaxation time spectrum model that reproduces the relaxation modulus measurements and is the best smoothed in the class of continuous square-integrable functions is sought. Based on the Hilbert projection theorem, the best-smoothed relaxation spectrum model is found to be described by a finite sum of specific exponential–hyperbolic basis functions. For noise-corrupted measurements, a quadratic with respect to the Lagrange multipliers term is introduced into the Lagrangian functional of the model’s best smoothing problem. As a result, a small model error of the relaxation modulus model is obtained, which increases the model’s robustness. The necessary and sufficient optimality conditions are derived whose unique solution yields a direct analytical formula of the best-smoothed relaxation spectrum model. The related model of the relaxation modulus is given. A computational identification algorithm using the singular value decomposition is presented, which can be easily implemented in any computing environment. The approximation error, model smoothness, noise robustness, and identifiability of the polymer real spectrum are studied analytically. It is demonstrated by numerical studies that the algorithm proposed can be successfully applied for the identification of one- and two-mode Gaussian-like relaxation spectra. The applicability of this approach to determining the Baumgaertel, Schausberger, and Winter spectrum is also examined, and it is shown that it is well approximated for higher frequencies and, in particular, in the neighborhood of the local maximum. However, the comparison of the asymptotic properties of the best-smoothed spectrum model and the BSW model a priori excludes a good approximation of the spectrum in the close neighborhood of zero-relaxation time.


Introduction
The relaxation spectrum is vital for constitutive models and for providing insight into the mechanical properties of polymers since, from the relaxation spectrum, other material functions used to describe rheological properties can be uniquely derived [1][2][3][4][5].It is applied for description, analysis, and to accomplish the pre-assumed mechanical properties of different polymers [3,[6][7][8].
The spectrum is not directly accessible via measurement and must be recovered from relaxation stress or oscillatory shear data.Numerous different methods have been proposed Polymers 2024, 16, 2300 2 of 33 during the last seven decades for relaxation spectrum identification using data both from the stress relaxation experiment [6,[9][10][11][12][13][14][15][16][17][18][19][20] and dynamic modulus tests [5,[21][22][23][24][25][26][27][28][29][30][31].The problem of relaxation spectrum identification is the ill-posed inverse problem of solving a system of Fredholm integral equations of the first kind obtained for discrete measurements of the relaxation modulus or the storage and loss modulus data.Therefore, the solutions, if any, are very sensitive to even small changes in the measurement data, which can lead to arbitrarily large changes in the determined relaxation spectrum.In consequence, robustly stable algorithms are required to solve it.The regularization of the original problem and/or constraining the set of admissible solutions is often necessary to construct such algorithms.
In the first works concerning the relaxation spectrum determination, the sets of the spectrum models were constrained to rather narrow classes of models.In 1948 Macey [9], while examining the viscoelastic properties of ceramic material, described the spectrum by the exponential-hyperbolic model, which corresponds to the modified Bessel function of the second kind and zero-order modeling the relaxation modulus.To describe the mechanical properties of polyisobutylene, Sips [10] introduced in 1950 a simple relaxation spectrum model given by the difference between two exponential functions, which implied a logarithmic model of the relaxation modulus.This model was next augmented to consider a long-term modulus by Yamamoto [11] and applied to study the rheological properties of the plant cell wall.The relaxation spectrum modeling in [9][10][11] is based on the known pairs of Laplace transforms.
The relaxation spectrum identification based on the Post-Widder formula [12] for the inverse Laplace transform was initiated by Alfrey and Doty [13], who proposed a simple differential model based on the first-order Post-Widder formula.Ter Haar [14] approximated the spectrum of relaxation frequencies using the modulus multiplied by time, the inverse of the relaxation frequency, which is, in fact, the Post-Widder inversion formula of the zero order.After many years, Bažant and Yunping [15] and Goangseup and Bažant [16] introduced a two-stage approach of approximating the stress relaxation data via multiple differentiable models of the relaxation modulus and, next, by applying the Post-Widder formula to designate the related model of the spectrum.The effectiveness of this approach depended, among other aspects, on the function applied to approximate the relaxation modulus.In [15], a logarithmic-exponential model of the relaxation modulus was proposed, for which the authors stated the third-order Post-Widder approximation to be satisfactory.
Both the algorithms based on the Post-Widder formula [13][14][15][16], using the least-squares approximation to guarantee the best fitting of the relaxation modulus measurement data, and those using the pairs of Laplace transforms [9][10][11], did not take into account the ill-posed nature of the relaxation spectrum determination problem.
Baumgaertel and Winter [21] used a nonlinear least-squares method for the recovery of a discrete relaxation time spectrum based on storage and loss modulus data, in which the number of discrete model modes was adjusted during the scheme iterations to avoid an illposedness of the problem and to enhance the model fit.Regularization was not applied here, as in several of the works discussed subsequently.Malkin [22] approximated a continuous relaxation spectrum using three constants: the maximum relaxation time, the slope in the logarithmic scale, and the form factor. Malkin et al. [23] derived a method of continuous relaxation spectrum calculations using the Mellin integral transform.The algorithm for the relaxation time spectrum approximation by power series was developed by Cho [24], which, using the regression of the dynamic modulus, provided a relaxation time spectrum as stable as the regularization method.The least squares identification without regularization was also applied by Babaei et al. [17] to determine the discrete Maxwell relaxation spectrum based on the stress relaxation data from the ramp test.Lv et al. [5] applied the extended least squares method (without regularization) to dynamic experiment data.Lee et al. [25] used the Chebyshev polynomials of the first kind to approximate dynamic moduli data and next derived a spectrum equation using the complex decomposition method and the Fuoss-Kirkwood relation without any regularization.Also, a derivative-based algorithm Polymers 2024, 16, 2300 3 of 33 for continuous spectrum recovery, which is also appropriate for the experimental situation where oscillatory shear data are only available for a finite range of frequencies, as proposed by Anderssen et al. [26], does not use regularization.
Honerkamp and Weese [27,28] combined nonlinear regression with Tikhonov regularization and proposed a specific viscoelastic model described by the two-mode log-normal function.In turn, Davies and Goulding [29] approximated the relaxation spectrum by a sum of scaling kernel functions located at appropriately chosen points.In the Mustapha and Phillips algorithm [30], the sequence of nonlinear regularized least-squares problems, solved with respect to both the discrete relaxation times and the elastic moduli, was performed with an increasing number of discrete model modes.The approach proposed by Stadler and Bailly [31] is based on the relaxation spectrum approximation using a piecewise cubic Hermite spline with respective regularization.The regularized algorithms presented in [27][28][29][30][31] were developed for dynamic rheological tests.A methodology to calculate the relaxation spectrum of biopolymeric materials from stress relaxation data has been proposed by Kontogiorgos [32], combining Hansen's least-squares numerical algorithm and Tikhonov regularization with the L-curve criterion chosen to select regularization parameter.Stankiewicz [18,19] and Stankiewicz et al. [20] derived different identification algorithms for the optimal regularized least-squares identification of relaxation time and frequency spectra in the classes of models defined by a finite series of different basis functions.
All the known methods for the relaxation spectrum optimal identification are based on the minimization of the quadratic model error defined directly for the measurements of the relaxation modulus or storage and loss modules.For example, in [5,17,21], the least-squares criterion was used directly, while in [18][19][20][21]27,28,30], regularized least-squares were used with various rules applied for the choice of regularization parameters to ensure the stability of the scheme and smoothness of the determined relaxation spectrum.In these papers, the mathematical formula describing the relaxation spectrum model was also, in advance, limited to the assumed class of admissible models.
In this paper, as a remedy for the ill-posed nature of the spectrum identification problem, the relaxation time spectrum model that reproduces the relaxation modulus measurements and which is the best smoothed in the class of continuous square-integrable functions was sought.This problem was formulated and solved in this paper for the first time for the spectrum of relaxation times.First, by applying the well-known Hilbert projection theorem, a new model was derived in which the best smoothing was achieved together with the simultaneous interpolation of relaxation modulus measurements.Next, to achieve noise robustness, the problem of the optimal smoothness of the spectrum model was augmented by introducing a quadratic term in the Lagrange functional of the original optimal spectrum smoothing problem.The necessary and sufficient optimality condition of the modified problem implied the best relaxation spectrum model as a finite sum of the basis functions given by the quotient of the exponential function and relaxation time.The components of the corresponding relaxation modulus model were given by simple hyperbolic functions.The permitted in advance small error of the relaxation modulus model combined with the specific modification of the Lagrange functional resulted in the model's noise robustness.The complete computational procedure for determining the best-smoothed model was given.The singular value decomposition method was used for algebraic computations.Analytical formulas describing the relaxation modulus model error, the relaxation spectrum smoothness, and noise robustness indices were derived as quadratic positive definite forms dependent on the sampling instants applied in the experiment and the relaxation modulus measurements.The monotonicity of these indices was analyzed.The applicability of the proposed model and algorithm to determining the optimally smoothed models of polymers characterized by the short and middle relaxation times of the Gauss-like relaxation spectra and the long relaxation times of the Baumgaertel, Schausberger, and Winter spectrum was verified.The rough applicability analysis of the proposed approach to modeling the relaxation time spectra of different types, such as the Kohlrausch-Williams-Watts, fractional Maxwell, Scott-Blair, inverse power, and multiplied power-exponential laws, was also carried based on the compatibility of the boundary conditions of the real spectra and the best-smoothed model.These studies have shown the applicability of the new model and identification algorithm for the optimal recovery of the smoothed relaxation spectrum of polymers with a very wide range of relaxation times.
In summary, this paper addresses the ill-posed problem of identifying the relaxation time spectrum in a new, original way, previously unknown in the literature.In the known methods, different mechanisms are built in to obtain a smooth enough and noise-robust relaxation spectrum model.Here, a new approach is proposed where the optimally smoothed continuous square-integrable model of the relaxation spectrum, which reproduces relaxation modulus measurements with assumed acceptable small errors, is directly sought.This problem is mathematically formulated and solved, resulting in a unique, best-smoothed relaxation spectrum model and a complete identification algorithm.In the construction of known relaxation spectrum identification methods, the primary idea was the best model approximation, and the next one, implied by the ill-posed nature of the task, was the concept of smoothing the model by regularizing the original problem of the model optimal approximation.Here, the main measure of the model's quality was the integral of the square of the relaxation spectrum, being simultaneously the measure of the model's smoothness.The idea of the optimal model approximation is replaced here by the classical interpolation of measurement points.In the basic problem, precise interpolation is applied and is next modified to interpolation with a small error being allowed to ensure the noise robustness of the model and algorithm.The idea of the Tikhonov regularization technique results from the essence of the ill-posed problem-the lack of uniqueness for its solution or its discontinuity with respect to the measurement data.The problem of smoothing the relaxation spectrum posed here finds inspiration in the consequences of the ill-posed problem and sometimes catastrophic fluctuations of the obtained solution, and it eliminates these model-devastating effects.
In Appendix A, the proofs of the main results and derivations of some mathematical formulas are given.Some tables of numerical results were moved to Appendix B to increase the clarity of the article.

Relaxation Time Spectrum
The continuous relaxation time spectrum H(τ) of a linear viscoelastic material is defined by the following integral [1,26]: where G(t) is the linear relaxation modulus at time t.The spectrum H(τ) is interpreted as a generalization of the discrete Maxwell spectrum to a continuous function [1,26] and characterizes the distributions of relaxation times τ.

Identification
A classical way of conducting the identification experiment studying viscoelasticity is the stress-relaxation test, where time-dependent stress is studied for the step increase in the strain [1,[34][35][36].Suppose a certain stress relaxation experiment resulted in a set of measurements of the relaxation modulus G(t i ) = G(t i ) + z(t i ) at the sampling instants t i > 0, i = 1, . . ., N, where z(t i ) is the measurement noise.Identification, classically, consists of the selection, within the chosen class of models, of such a model that ensures the best fit to the measurement results.As a measure of the model's accuracy, the quadratic index, used in the least squares approach, is usually taken.However, in this paper, we look, in the class of continuous square-integrable functions, for the best-smoothed model H M (τ) that reproduces relaxation modulus measurements G(t i ).This problem is formulated and solved in the next section.

Results and Discussion
In this section, the problem of optimal smoothing of the relaxation spectrum model is mathematically formulated and solved using the Hilbert projection theorem.As a result, the best-smoothed relaxation spectrum model is derived in the form of a series of specific basis functions given by the quotient of the exponential function and the relaxation time.The respective model of the relaxation modulus was found to be described by a series of hyperbolic functions.The properties of these basis functions were examined, and the identifiable property of the best-smoothed model was demonstrated.For noise measurement data, the modification of the problem for the spectrum model's smoothness was proposed by augmenting the Lagrange functional.A dual approach was applied to solve the modified problem, resulting in the necessary and sufficient optimality condition for the optimal relaxation spectrum model.Direct analytical formulas for the relaxation spectrum and modulus models are given; their numerical realization by the singular value decomposition of the basic matrix was proposed.The model smoothness, noise robustness for noisy measurements of the relaxation modulus, and the error of the relaxation modulus model were analyzed.A simple identification scheme was proposed.Finally, the results of simulation studies for polymers described by Gaussian-like and BSW relaxation spectra distributions are presented.

The Problem of Optimal Smoothing of the Relaxation Spectrum Model
Consider the following problem.Find function H M (τ) ∈ L 2 (0, ∞) that minimizes the integral square index: under the constraints Note, that the set of functions H M (τ) ∈ L 2 (0, ∞) satisfying linear constraints (3) is closed and convex.Since the Hilbert projection theorem [33] implies the existence of a unique element with a minimal norm in the nonempty closed and convex subset of the Hilbert space, the existence and uniqueness of the solution to the smoothing problem (2) and (3) follows.
The Lagrange functional of the optimization problem defined by (2) and ( 3) is as follows: where λ N = [λ 1 , . . . ,λ N ] T is a vector of Lagrange multipliers λ i .The necessary and suffi- cient optimality conditions for the linear-quadratic optimization task (2) and (3) are given by the equation Polymers 2024, 16, 2300 6 of 33 together with the constraints of (3).Substituting H M (τ), given by (5), into (3) yields the following system of equations: By applying the substitution v = 1/τ in the integrals of the right-hand side of (6) we obtain: Introducing the N × N dimensional symmetric matrix composed of the elements Φ ij in the i row and j column according to the system of Equation ( 6) can be rewritten in compact form as with the vector of the relaxation modulus measurements where superscript "T" indicates transpose.The main properties of the matrix Φ N are summarized in the following proposition shown in Appendix A.1.
Proposition 1.For arbitrary sampling instants t i > 0, i = 1, . . ., N, such that t i+1 > t i , a symmetric N × N matrix Φ N defined by (7) and ( 8) is a positive definite Gram matrix, which can be expressed as is the unique symmetric non-singular positive definite square root of Φ N .Then, the inverse matrix N is a positive definite too.
Since matrix Φ N is non-singular, the unique solution of ( 9) is given by Therefore, by virtue of (5), the best-smoothed model is described by the finite series of the basis functions (compare (5)) The lower index 'N' is used in H N (τ) to express the dependence on the number of measurements.
Introducing the notation Polymers 2024, 16, 2300 and bearing in mind (11), model (12) can be expressed in compact vector-matrix form as In view of (1) and ( 12), the related model of the relaxation modulus is as follows: and can be described by the finite series of hyperbolic basis functions (compare (7)) Similarly, using the right equality in (16) and Equation ( 11), we obtain where, in view of ( 17), the vector function φ N (t) is as follows: The index minimized in (2), is a direct measure of smoothing the relaxation spectrum model.For the optimal model H N (τ) (15), the smoothness index is as follows: Using ( 14), (7), and ( 8) we find Thus, Equation ( 20) yields which means that the smoothness of the optimal model depends both on the time instants t i selected for the stress relaxation experiment, affecting the matrix Φ N and on the experiment results G N .As a result, the following result can be stated.
Theorem 1.For arbitrary sampling instants t i > 0, i = 1, . . ., N, such that t i+1 > t i , the unique optimally smoothed model of the relaxation time spectrum defined by the optimization task (2) and ( 3) is given by H where the vector functions h N (τ) and φ N (t) are defined by ( 14) and ( 18), respectively, and matrix Φ N is defined by ( 8) and (7).
The relaxation spectrum H N (τ) that solves the optimization task ( 2) and (3), is the most smoothed model in the class of square-integrable functions that, simultaneously, guarantee the reconstruction of the measurements of the relaxation modulus.Some useful algebraic identities concerning the matrix Φ N and vector φ N (t) are given in Appendix A.2.

Properties of the Basis Functions
The basis functions h i (τ) (13) and φ i (t) (17) of the relaxation spectrum and modulus models are positive definite and depend on the times t i applied in the stress relaxation experiment.The greater the sampling instants t i , the faster the basic functions φ i (t) decrease.By (13), the first derivative is as follows: Thus, the basis functions h i (τ) for τ = t i have a global maximum equal to h i (τ) = 1/(et i ), which decreases with increasing index i due to the assumed monotonicity of the sequence {t i }.This means that increasing the number of measurements N, i.e., increasing the model components, can allow for the good modeling of multimodal spectra, which is confirmed by the second example presented in the final part of this paper.
Since for τ → 0 + , using the L'Hospital's rule, we have and for τ → ∞ the functions are h i (τ) → 0 , the best model H N (τ) tends to zero both for τ → 0 + and τ → ∞ (zero boundary conditions).The basis functions φ i (t) given by (17)  monotonically decrease to zero as t → ∞ .
The five first basis functions h i (τ) (13) are shown in Figure 1 for the sampling instants t i = 10, 30, 50, 70, 90 and t i = 0.1, 0.5, 1, 1.5, 2 s. Figure 2 shows the related functions φ i (t) (17).The logarithmic scale is applied for the time axes in these figures.The basis functions h i (τ) and φ i (t) are expressed in s −1 .From Figure 2, it can be seen that the monotonicity of basis functions φ i (t) is in good agreement with the courses of the relaxation modulus obtained in an experiment for real polymers; for example, these include elastic polyacrylamide hydrogels [35] (Figures 2a,b, 4a, A5, A7 and A8a), concrete [37] (Figure 13) and rubber [38] (Figure 2).ation experiment.The greater the sampling instants  , the faster the basic functions   decrease.By (13), the first derivative is as follows: Thus, the basis functions ℎ  for   have a global maximum equal to ℎ  1  ⁄ , which decreases with increasing index  due to the assumed monotonicity of the sequence  .This means that increasing the number of measurements , i.e., increasing the model components, can allow for the good modeling of multimodal spectra, which is confirmed by the second example presented in the final part of this paper.
Since for  → 0 , using the L'Hospital's rule, we have and for  → ∞ the functions are ℎ  → 0, the best model   tends to zero both for  → 0 and  → ∞ (zero boundary conditions).The basis functions   given by ( 17) monotonically decrease to zero as  → ∞.
The five first basis functions ℎ  (13) are shown in Figure 1 for the sampling instants  10,30,50,70,90 and  0.1, 0.5, 1, 1.5, 2 s. Figure 2 shows the related functions   (17).The logarithmic scale is applied for the time axes in these figures.The basis functions ℎ  and   are expressed in  .From Figure 2, it can be seen that the monotonicity of basis functions   is in good agreement with the courses of the relaxation modulus obtained in an experiment for real polymers; for example, these include elastic polyacrylamide hydrogels [35] (Figures 2a,b, 4a, A5, A7, and A8a), concrete [37] (Figure 13) and rubber [38] (Figure 2).

Identifiability
The basic and obvious requirement for any identification method is that if the real characteristic is described by a model from the considered class of models, and the measurements are noise-free, then the method should guarantee the unique determination of the real characteristic, i.e., ensure its identifiability [39,40].
Assume that the real spectrum is of the form where  represents the real parameters.Introducing the vector   , … ,  and bearing in mind ( 13) and ( 14), spectrum (23) can be expressed as Assume that the measurements of the relaxation modulus are noise-free.Thus, for  1, … , , by virtue of (1), ( 7) and ( 23), we have Therefore, using ( 8) and (10), the vector of the relaxation modulus measurements can be expressed as whence, according to ( 15) and ( 24), the relaxation spectrum model is as follows: i.e., the model smoothing identification results in the determination of the real relaxation spectrum (23).

Modification
The value of the Lagrange multiplier  is the dual price [41], which in problem (2) and ( 3) is "paid" for satisfying the -th constraint,  1, … , .The higher the value of  (precisely, the modulus of  ), the more difficult it is to meet this constraint and the stronger the chains it imposes.The impact of the fluctuations in the measurements of the relaxation modulus  ̅  , i.e., the impact of changes in the left side of each of Equation (3), on the smoothness of the spectrum is then greater.Really, by (22) and ( 11), we have Basis functions φ i (t) (17), i = 1, . . ., 5, of the relaxation modulus model G N (t) (16) for time sampling instants: (a) t i = 10, 30, 50, 70, 90 s and (b) t i = 0.1, 0.5, 1, 1.5, 2 s.

Identifiability
The basic and obvious requirement for any identification method is that if the real characteristic is described by a model from the considered class of models, and the measurements are noise-free, then the method should guarantee the unique determination of the real characteristic, i.e., ensure its identifiability [39,40].
Assume that the real spectrum is of the form where a j represents the real parameters.Introducing the vector a = [a 1 , . . . ,a N ] T and bearing in mind ( 13) and ( 14), spectrum (23) can be expressed as Assume that the measurements of the relaxation modulus are noise-free.Thus, for i = 1, . . ., N, by virtue of (1), ( 7) and ( 23), we have Therefore, using ( 8) and (10), the vector of the relaxation modulus measurements can be expressed as whence, according to ( 15) and ( 24), the relaxation spectrum model is as follows: i.e., the model smoothing identification results in the determination of the real relaxation spectrum (23).

Modification
The value of the Lagrange multiplier λ i is the dual price [41], which in problem (2) and ( 3) is "paid" for satisfying the i-th constraint, i = 1, . . ., N. The higher the value of λ i (precisely, the modulus of λ i ), the more difficult it is to meet this constraint and the stronger the chains it imposes.The impact of the fluctuations in the measurements of the relaxation modulus G(t i ), i.e., the impact of changes in the left side of each of Equation (3), on the smoothness of the spectrum is then greater.Really, by ( 22) and ( 11), we have Therefore, the vector of the optimal Lagrange multipliers is the measure of the model's smoothness index sensitivity with respect to the fluctuations of the relaxation modulus measurements.To reduce it, the value of the multiplier λ i should be reduced; precisely, the values of the modulus of λ i should be reduced.
The Lagrange functional of the optimization problem defined by ( 2) and ( 3) is described in Equation ( 4).In order to decrease the values of |λ i | and the functional (4), being maximized with respect to λ i according to the dual approach, is modified by introducing the quadratic term ∑ N i=1 λ 2 i = λ T N λ N , which means that the modified Lagrange functional is defined as follows: where γ is a small positive constant and the weight that represents the relative importance of the square component λ T N λ N with respect to the original (non-modified) Lagrange functional L(H M , λ N ) (4).The parameter γ has no physical interpretation as the regularization parameter in the classical Tikhonov regularization.
In Appendix A.3, the unique saddle point of the modified Lagrange functional ( 25) is found.The saddle point defines the modified best-smoothed model H γ N (τ), which depends on the measurements t i , G(t i ) and the parameter γ introduced in (25).Theorem 2. For arbitrary sampling instants t i > 0, i = 1, . . ., N, such that t i+1 > t i and the arbitrary non-negative parameter γ, the model of the relaxation time spectrum H γ N (τ) defined by the unique saddle point of the modified Lagrange functional (25) is given by while the respective relaxation modulus model where the vector functions h N (τ) and φ N (t) are defined by ( 14) and ( 18), respectively, matrix Φ N is defined by ( 7) and ( 8) and I N,N is the N × N dimensional unit matrix.
The upper index γ in the notations indicates the dependence on the parameter γ introduced in the modified Lagrange functional (25).
To achieve the dimensional homogeneity of the components of the Lagrange functional (25), the multipliers λ i are expressed in Pa•s, while the unit of the parameter γ is s −1 .The dimensional homogeneity of the matrix Φ N + 4γI N,N is then achieved.
It is demonstrated in Appendix A.4 that for the optimal vectors of the Lagrange multipliers, λ N in (11) for original optimization task (2) and (4) and the vector λ γ N (A13) of the saddle point of the modified Lagrange functional (25), the following inequality holds: which means that λ γ N < ∥λ N ∥, i.e., the purpose of the modification introduced into the Lagrange functional at the beginning of this section was achieved.

Model Error
Let us introduce, by analogy to the measurements vector G N (10), the vector of the values of the relaxation modulus model ( 27) for all sampling instants t i : For any t = t i , by virtue of (27), whence, bearing in mind Equation (A4), we have Therefore, for the model parameterized by γ > 0, the relaxation modulus equations in (3) are not satisfied and the error of these equations is as follows: Through the identity (A6), the model error ε N can be expressed as and, bearing in mind (A13), can be equivalently expressed as ε N = 2γλ γ N .Therefore, for γ = 0 the model error ε N = 0 N , which is clear since in the original task ( 2) and (3) the constraints in Equation (3) are exactly satisfied; here, 0 N denotes the N dimensional vector of zero elements.
By (30), the square model error is given as follows: In Appendix A.5, the following result is proved.
Proposition 2. For arbitrary sampling instants t i > 0, i = 1, . . ., N, such that t i+1 > t i , and the arbitrary non-negative parameter γ, the square model error ε T N ε N (31) of the relaxation modulus equations monotonically increases as the function of the parameter γ , which is strictly convex for γ such that Φ N ≥ 8γI N,N , and strictly concave in the case Φ N ≤ 8γI N,N .

Smoothness
For model H γ N (τ) (26), the smoothness index I N (19) is as follows: and, bearing in mind (21), can be expressed as or, by applying identity (A7), an equivalent form useful for further differential analysis can be obtained Polymers 2024, 16,2300 Since, according to (22), the smoothness index for the original optimization task (2) and ( 3) can be expresses as N G N , by inequality (A5), the next estimation follows: I γ N < I N for any γ > 0 and arbitrary measurement data, i.e., the smoothness of the spectrum model H γ N (τ) (26) is stronger that the smoothness of the original model H N (τ) (15); this was the idea of the modification introduced in the Lagrange functional.
In Appendix A.6, the following formulas describing the first and second derivatives of I γ N with respect to γ are derived: and Thus, the smoothness index I γ N is the monotonically decreasing convex function of the parameter γ > 0. The following rule holds: the greater the parameter γ is, the more highly bounded the fluctuations of the spectrum model H γ N (τ) (26) are.

Noise Robustness
Following [18,19], as a reference point for the model H γ N (τ) described by Equation (26), the model of the spectrum that can obtain for the same parameter γ, the same number of measurements N and the same time instants t i on the basis of ideal measurements of the relaxation modulus is considered, which is described as follows: where G N is the vector of the noise-free relaxation modulus, i.e., In view of ( 26) and (36), we find where T is the vector of measurement noises.
Consider the square integral index By (37), this index can be expressed as whence, by virtue of (21), we obtain which, using the Gram property of Φ N and using the identity (A7), can by rewritten as follows: Therefore, the noise robustness depends on the parameter γ, the measurement noises and the sampling instants t i that uniquely determine the matrix Φ N .Since both models are continuous with respect to the relaxation time τ, by virtue of (40), for any non-negative γ, the spectrum H γ N (τ) tends to the noise-free spectrum ∼ H γ N (τ) for each time τ linearly with respect to the norm ∥z N ∥, as ∥z N ∥ → 0 , and the faster this is, the larger the parameter γ.
By (A5), we have which means better noise robustness than for the original model H N (τ) (15) for any γ > 0.
Similarly, as for the smoothness index I γ N (compare indices (33) and ( 40)), the following formulas describing derivatives of Q γ N with respect to γ were derived: and which means that Q γ N is the monotonically decreasing convex function of the parameter γ > 0 that takes the maximal value equal to Q 0 N (41) for γ = 0.

Algebraic Background of the Computational Algorithm
The singular value decomposition (SVD, [42]) technique can be used in numerical computations in order to determine the inverse matrix (Φ N + 4γI N,N ) −1 in (26).Let SVD of the N × N dimensional matrix Φ N (8) take the following form [42]: where Σ N = diag(σ 1 , . . ., σ N )ϵR N,N is the diagonal matrix containing the singular values σ i of the matrix Φ N [42], and the matrix U N ∈ R N,N is orthogonal.Thus, where the N × N diagonal matrix (Σ N + 4γI N,N ) −1 is as follows: Therefore, the optimal spectrum model ( 26) can be described by while the respective model (27) of the relaxation modulus is as follows: where the vector of model parameters is Using ( 42) and ( 43), the smoothness index I γ N (32) is expressed as where the N × N diagonal matrix Similarly, using ( 42) and ( 43), the noise robustness index Q γ N (39) can be rewritten as Combining (31) and (43), we obtain the formula describing the square model error

Algorithm
The determination of the best-smoothed model of the relaxation time spectrum involves the following steps: 1.
Compute the vector of model parameters g γ N (45).
Determine the square model error ε T N ε N according to (47) and the smoothness index I γ N using Formula (46).7.
Check if the smoothness of the spectrum model H Only the SVD of the matrix Φ N of computational complexity O N 3 [42] is a spaceand time-consuming task in the scheme.However, for given sampling points, the SVD must be computed only once in step 3.The matrix Φ N does not depend on the relaxation modulus measurements G(t i ).Therefore, when the identification scheme is applied for successive samples of the same material, step 3 should not be repeated whenever the same time instants t i are kept in the experiment.This is because, using (44) and (45), we have where the vector function depends only on the sampling points t i and does not depend on the relaxation modulus measurements G(t i ).Therefore, the function ϑ γ N (τ) must be computed only once and used to determine the model H γ N (τ) for many samples whenever the same t i is kept.
The best-smoothed spectra models, the values of the square model error ε T N ε N (31), where the error ε N is defined by (29), the smoothness index given by Formula (32), and the square noise robustness index Q γ N (38) expressed by Equation ( 39) are examined for different number of the measurements N and different values of the parameter γ.
The "real" materials and the best-smoothed models were simulated in Matlab R2023b, The Mathworks, Inc., Natick, MA, USA.For the singular value decomposition procedure svd was applied.A normal distribution with zero mean value and variance σ 2 as well as the uniform distribution were applied for the random independent generation of additive measurement noises.

Example I
Considering the polymer whose relaxation spectrum is described by the uni-modal Gaussian-like distribution as follows: and where the parameters are as follows [54,55]: ϑ = 31, 520 Pa•s, m = 0.0912 s −1 and q = 3.25 × 10 −3 s −2 .The related relaxation modulus is desribed by the function [55]: where the complementary error function er f c(x) is given by [56]: In the experiment, N sampling instants t i were generated with the constant period in the time interval of T = [0, 200] seconds, selected on the basis of the modulus G(t) (49) course.

Noise-Free Measurements
For noise-free measurements of the modulus G(t) the best-smoothed model solving the original task (2) and (3) was determined for N = 20, 100, 150, 200, 500, 1000 measurements.Two optimal models H N (τ) (15) and the 'real' spectrum H(τ) (48) are plotted in Figure 3. Small subfigures confirm the excellent model fit; real spectra described by the red lines practically coincide with the blue models both for the small and large number of measurements.This shows that in the case of noise-free measurements, the practically ideal approximation of the real relaxation spectrum was obtained even for a small number of measurements (N = 20).In Figure 4, the related models of the relaxation modulus G N (t) (16) are plotted; the measurements G(t i ) of the 'real' modulus G(t) (49) are also marked.The values of the smoothness index I N (19) are given in Table 1.red lines practically coincide with the blue models both for the small and large number of measurements.This shows that in the case of noise-free measurements, the practically ideal approximation of the real relaxation spectrum was obtained even for a small number of measurements ( 20).In Figure 4, the related models of the relaxation modulus  ̅  (16) are plotted; the measurements  ̅  of the 'real' modulus   (49) are also marked.The values of the smoothness index ℐ (19) are given in Table 1.
(a) (b)     Additive independent measurement noises are generated by a normal distribution with zero mean value and variance  .For the noise robustness analysis, the standard deviations  2, 4, 6 Pa were used.The parameters  5 10 , 10 , 5 10 , 10 s were applied.
In Table 2, the values of the square model error   (31), the smoothness index ℐ (32), and the square noise robustness index  (38) are given for noises of  2 Pa, while for the stronger noises, the same data are given in Tables A1 and A2 in Appendix B. As previously, the exemplary courses of the spectrum models   (26) for  20 and  500 measurements are illustrated in Figure 5 for noises of  2, 4, 6 Pa, while the respective relaxation modulus models  ̅  (27) are depicted in Figure 6.
An inspection of Figure 5a,c,e shows that for each noise case, the number of  20 measurements was not enough to obtain the satisfactory smoothness of the model   even for the weakest noises.So,  20, which is good for noise-free case, fails here.However, for  500 measurements, the model   is smoothed enough, and the
In Table 2, the values of the square model error ε T N ε N (31), the smoothness index I γ N (32), and the square noise robustness index Q γ N (38) are given for noises of σ = 2Pa, while for the stronger noises, the same data are given in Tables A1 and A2 in Appendix B. As previously, the exemplary courses of the spectrum models H γ N (τ) (26) for N = 20 and N = 500 measurements are illustrated in Figure 5 for noises of σ = 2, 4, 6 Pa, while the respective relaxation modulus models G γ N (t) (27) are depicted in Figure 6.

Table 2.
For Example I, the square model error ε T N ε N (31), the smoothness index dτ, Equation ( 32), and the noise robustness index Q γ N (38) for N measurements of the relaxation modulus corrupted by normally distributed additive independent noises with zero mean value and standard deviation σ = 2 Pa; parameter γ is introduced in the modified Lagrange Functional (25).An inspection of Figure 5a,c,e shows that for each noise case, the number of N = 20 measurements was not enough to obtain the satisfactory smoothness of the model H γ N (τ) even for the weakest noises.So, N = 20, which is good for noise-free case, fails here.However, for N = 500 measurements, the model H γ N (τ) is smoothed enough, and the influence of the regularization parameter γ is much weaker; see Figure 5b,d,f.Figure 6 and the values of the model errors ε T N ε N from Tables 2, A1 and A2 confirm the excellent approximation of the relaxation modulus model, even though model imbalance is allowed.
The analytically shown monotonicity of the smoothness I γ N and noise robustness Q γ N indices, being the monotonically decreasing convex functions of the parameter γ, is reflected in the numerical studies.An inspection of the numerical results indicates that for N ≥ 50 and any fixed parameter γ, the smoothness index I γ N is a monotonically increasing function of the number of measurements; the slower this is, the larger the number N.An analysis of the asymptotic properties of the algorithm and optimal model H γ N (τ) (26) as the number of measurements grows to infinity will be the subject of future studies.

Example II
Consider the double-mode Gaussian-like distribution of the relaxation spectrum [19,20,44] where the parameters are as follows: The double-Gaussian relaxation spectra are examined while developing new identification methods in [31] (Figure 2), ).Such spectra describe the rheological properties of various polymers [44] (Figures 4b and 8b), polyacrylamide gels [35] (Figure A4), and wood [38].The corresponding 'real' relaxation modulus is composed of two summands described by formulas like that of (49).In the experiment, N = 50, 100, 200, 500, 1000, 5000 sampling instants t i were generated with the constant period in the time interval T = [0, 1550] seconds, selected in view of the course of the mod- ulus.Following [19,20], additive measurement noises z(t i ) were selected independently by random choice with uniform distribution on the interval [−0.005, 0.005] Pa.In Table 3, the values of the square model error ε T N ε N (31), the smoothness index I γ N (32), and the square noise robustness index Q γ N (38) are given.The spectrum models H γ N (τ) (26) are illustrated in Figure 7 along with the real spectrum (50).Since, similar to the one-mode Gaussian relaxation spectrum, the relaxation modulus models G γ N (t) (27) for different N and γ values practically coincide, the respective figures are omitted here.Table 3.For the relaxation spectrum (50) from Example II described by double-mode Gaussian distribution: the square model error ε T N ε N (31), the smoothness index I γ N (32), and the noise robustness index Q γ N (38) for N measurements of the relaxation modulus corrupted by additive independent noises uniformly distributed on the interval [−0.005, 0.005] Pa; parameter γ introduced in the modi- fied Lagrange functional (25).
Additive measurement noises z(t i ) were selected independently by random choice with uniform distribution on the interval [−0.005, 0.005] MPa.The results of the numerical experiment are given in Table 4 and illustrated in Figure 8.  (38) for N measurements of the relaxation modulus corrupted by additive independent noises selected according to uniform distribution from the interval [−0.005, 0.005] MPa and parameters γ introduced in the modified Lagrange functional (25).Since the real spectrum H(τ) (51) tends to infinity for τ → 0 whenever at least one of the parameters ρ 1 and ρ 2 is negative, the best-smoothed model H γ N (τ) (26) cannot adequately approximate this spectrum for small relaxation times τ; in the example for 0 < τ < 10 3 s.This is well illustrated in Figure 8.However, this figure also shows that for a sufficiently large γ, the spectrum H(τ) for higher frequencies and its local maximum are well approximated.

Applicability for Identification of Relaxation Spectra of Different Types
The natural condition of this approach's successful applicability follows from the properties of the best-smoothed model   (26) yielded by the properties of the basis functions ℎ  (13), which compose the vector   according to (14).Since for  → 0 and  → ∞, the basis functions are ℎ  → 0 (c.f., Section 3.2), the best model   also tends to zero as the relaxation time  tends to zero and to infinity.Therefore, zero boundary conditions limit the scope of applicability of the model and method to the real relaxation time spectra that satisfy these conditions.The example of the BSW spectrum demonstrates that the properties of the spectrum for  → 0 are essential here since the real relaxation time spectra and the known spectra models tend to zero as the relaxation time  tends to infinity.

Applicability for Identification of Relaxation Spectra of Different Types
The natural condition of this approach's successful applicability follows from the properties of the best-smoothed model H γ N (τ) (26) yielded by the properties of the basis functions h i (τ) (13), which compose the vector h N (τ) according to (14).Since for τ → 0 + and τ → ∞ , the basis functions are h i (τ) → 0 (c.f., Section 3.2), the best model H γ N (τ) also tends to zero as the relaxation time τ tends to zero and to infinity.Therefore, zero boundary conditions limit the scope of applicability of the model and method to the real relaxation time spectra that satisfy these conditions.The example of the BSW spectrum demonstrates that the properties of the spectrum for τ → 0 + are essential here since the real relaxation time spectra and the known spectra models tend to zero as the relaxation time τ tends to infinity.
The Kohlrausch-Williams-Watts (KWW) model of the stretched exponential relaxation described by [57] where the stretching exponential 0 < β < 1, τ KWW is the characteristic relaxation time and G 0 denotes the initial shear modulus, has been found by many researchers to be more appropriate than standard exponentials [57][58][59][60][61][62].In spite of the simple, compact form of ( 52), the related unimodal [58] relaxation spectrum is described by the following infinite series [57,58]: which is based on Pollard's representation of the stretched exponential as a Laplace integral [63], where Γ(n) is Euler's gamma function [64] (Equation (A.1.1)).However, for some specific stretching exponentials, namely β = 1 2 , β = 1 3 and β = 2 3 , the KWW spectrum has a compact form described by some special functions [58].For (53), both zero boundary conditions are satisfied; compare [57] (Figure 1a).Therefore, the proposed approach can be used to identify the spectrum of materials whose relaxation processes are described by the KWW model, e.g., polymer melts [59], the local segmental dynamics of poly(vinylacetate) [60], the segmental dynamics and the glass transition behavior of poly(2-vinylpyridine) [61], the relaxation of bone and bone collagen [62], alginate films while considering glycerol concentration [65], and even the relaxation processes of the onion structure in sine-oscillatory shear [66].The best-smoothed model H γ N (τ) (26) given by finite series may prove to be more useful than the original KWW spectrum (53).
In recent decades, non-integer order differential equations have increased interest in the modeling of time-dependent relaxation processes; the fractional Maxwell model (FMM) and the elementary Scott-Blair model are probably the most known rheological non-integer order models.The applicability of the FMM relaxation time spectrum, which is described by the compact analytical formula [54] (Equation ( 12)), to modeling the unimodal relaxation spectra of polymers was recently examined in [54].However, it was demonstrated in [54] that the FMM relaxation spectrum tends to infinity as τ → 0 + [54] (Proposition 2, Equation ( 19)); therefore, the exact fitting of the FMM-type spectrum by the proposed best-smoothed model in the whole relaxation time domain is excluded, which is similar for the BSW spectrum.The relaxation time spectrum of the Scott-Blair model described by the inverse power of the relaxation time with the non-integer exponent, see [54] (Equation ( 15)), also loses the zero boundary condition at zero relaxation time.
Similarly, real relaxation spectra which are well characterized by simple inverse power laws with various exponents [67]; for example, the power-type spectrum of cross-linking polymers at their gel point described by Winter and Chambon model with an exponent of −1/2 [68] and the spectrum of solution-polymerized styrene butadiene rubber described by a combined four-interval power model with fractional exponents [69] could not be successfully identified by the proposed method in the whole relaxation times domain.In turn, Winter's power law relaxation time spectrum with a positive exponent [70,71] (Equation (2)), which was proposed to describe relaxation in many molecular and colloidal glasses, although satisfies the zero initial condition, could not be determined by the proposed algorithm due to its confined domain.
Polymers 2024, 16,2300 However, the proposed approach can be successfully applied to identify the relaxation spectra of materials such as bitumen, being characterized by the broadened power law model [71] (Equation ( 8)): which multiplicative form combines the power law with an exponential of stretching parameter β.In the above model, the exponent 0 < n α < 1, τ α is the longest relaxation time and G c is the plateau modulus.The unimodal spectrum (54) satisfies both zero boundary conditions; compare [71] (Figure 11a).

Conclusions
The objective of this paper was to develop a relaxation time spectrum model that could reproduce the relaxation modulus measurements and which is the best-smoothed in the class of continuous square-integrable functions.The unique optimal relaxation spectrum model was found to be described by a finite series of specific exponential-hyperbolic functions.A new identification algorithm was proposed in which the best smoothing of the model was achieved together with the simultaneous reconstruction of relaxation modulus measurements with small model errors.The analytical and numerical studies proved that using a developed model and algorithm, it is possible to determine the relaxation spectrum model for a wide class of polymers with zero boundary conditions, in particular, Gaussian-like distributed relaxation spectra.The model is smoothed and noise robust; small relaxation modulus model errors are obtained.The applicability of this approach to determining the Baumgaertel, Schausberger, and Winter spectrum was also examined, and it was proved that, due to the asymptotic properties of this spectrum, it can be well approximated for higher frequencies and, in particular, in the neighborhood of the local maximum.The rough applicability analysis, based on the consistency of the zero boundary conditions of the real spectra and the best-smoothed model, shows the possibility of using the proposed method and model to describe the relaxation spectra of different types that are characteristic of many polymers.However, the search for such a modification of the proposed approach so that it can also be applicable to the identification of spectra with non-zero boundary conditions for relaxation times approaching zero, like the BSW spectrum, will be the subject of future research.Generally, the properties of the method, including the smoothness of the relaxation spectrum, depend both on the experiment plan, i.e., on the sampling instants used in the relaxation test, and on the relaxation modulus measurements.The results of numerical studies confirm the analytically proved monotonic dependence on the gamma parameter: monotonically increasing for the square model error and monotonically decreasing for the noise robustness and spectrum model smoothness indices.However, the dependence of these indices on the number of measurements is not so clear and must be the subject of further studies.
Summarizing the numerical studies implies the following directions for future research: • The asymptotic analysis of the model and identification algorithm properties as the number of measurements tends to infinity; • The modification of the proposed approach for smoothing the spectrum model with a non-zero boundary condition for zero relaxation time; • The modification for non-zero equilibrium modulus; • The recurrent realization of the algorithm.This method can be applied for any deformation process described by definitional Equation (1), i.e., both for uniaxial deformation, uniaxial stress, and uniaxial stretching, assuming that the relaxation modulus of the respective process is experimentally accessible.The relaxation time spectrum in the respective state (uniaxial deformation, stress, or stretching) is then determined.An appropriate modification of the algorithm can be developed to apply the concept of optimal relaxation spectrum model smoothing for oscillatory shear measurements of the storage and loss moduli.This will be the subject of further research.
dτ, Equation (32), and the noise robustness index Q γ N (38) for N measurements of the relaxation modulus corrupted by normally distributed additive independent noises with zero mean value and the standard deviation σ = 4[Pa]; parameter γ is introduced in the modified Lagrange Functional (25).

γN
(τ) measured by I γ N and the error of the relaxation modulus model G γ N (t) measured by ε T N ε N are, simultaneously, satisfactory.If not, increase the parameter γ and repeat the computations starting from step 4. If yes, accept the current H γ N (τ) as the best-smoothed relaxation spec- trum model.

Figure 4 .
Figure 4.The measurements G(t i ) of the 'real' relaxation modulus G(t) (49) (red circles) from Example I and the optimal model G N (t)(16) for N noise-free relaxation modulus measurements: (a) N = 20; (b) N = 1000.

Figure 6 .Figure 5 .Figure 5 .
Figure 6.The measurements  ̅  of the 'real' relaxation modulus   (49) (red circles) from Example I and the model  ̅  (27) for  measurements of the relaxation modulus corrupted by

Figure 6 .
Figure 6.The measurements  ̅  of the 'real' relaxation modulus   (49) (red circles) from Example I and the model  ̅  (27) for  measurements of the relaxation modulus corrupted by

Figure 7 .
Figure 7. Relaxation time spectrum   (50) (solid red line) from Example II and the corresponding models   (26) determined for  measurements of the relaxation modulus corrupted by additive independent noises uniformly distributed on the interval 0.005, 0.005 Pa: (a)  50; (b)  100; (c)  200; (d)  500; (e)  1000; and (f)  5000; the values of the parameter  introduced in the modified Lagrange functional (25) are given in the figures.

Figure 8 .
Figure 8.The relaxation time BSW spectrum (51) (solid red line) from Example III and the corresponding models   (26) for  measurements of the relaxation modulus corrupted by additive independent noises selected according to uniform distribution from the interval 0.005, 0.005 MPa and parameters  introduced in the modified Lagrange functional (25): (a)  50; (b)  100; (c)  200; (d)  500; (e)  1000; and (f)  5000.

2 dτ,
Equation(32), and the noise robustness index Q γ N(38) for N measurements of the relaxation modulus corrupted by normally distributed additive independent noises with zero mean value and the standard deviation σ = 6[Pa]; parameter γ is introduced in the modified Lagrange Functional(25).

Table 1 .
(48)smoothness index I N(19)for noise-free N measurements of the relaxation modulus from Example I described by the relaxation spectrum H(τ)(48).

Table 1 .
The (48)thness index ℐ (19) for noise-free  measurements of the relaxation modulus from Example I described by the relaxation spectrum  (48).

Table 2 .
(25)Example I, the square model error  (31), the smoothness index ℐ   , Equation(32), and the noise robustness index (38)for  measurements of the relaxation modulus corrupted by normally distributed additive independent noises with zero mean value and standard deviation  2 Pa; parameter  is introduced in the modified Lagrange functional(25).